C
C     展开dipole field的贡献，理解如何展开
C     那么便知晓如何考虑小尺度的贡献
C
C  ------------------------------------------------
C  Calculate the dipole field at the positions of each sphere,
c  The pp0 and qq0 here are the corresponding quantities from xuy95
c  p0 and q0 are
c  multiplied by the Mie coefficients an and bn respectively,
c  compare xuy95 eqs. 14, 15, 30, 33, 85. Note that a different
c  normalization Emn is used here and in xuy95; see the gmm01f manual
c  for details.
c  ------------------------------------------------
c  set all as, bs, as0, bs0, p0, q0 to zero
C
C     模拟dipole作为入射场，VSH展开成平面波近似
C
      subroutine incident_dipole(nL, p0, q0, as, bs, pp0, qq0, aMie,
     $     bMie, atr0, btr0, ek, drot, nmax, uvmax, k, r0, fint)
      implicit double precision (a-h,o-z)
      include 'gmm01f.par'
      parameter (nmp=np*(np+2),nmp0=(np+1)*(np+4)/2)
      parameter (ni0=np*(np+1)*(2*np+1)/3+np*np)
      parameter (nij=nLp*(nLp-1)/2)

      integer nmax(nLp), uvmax(nLp), ind(nLp)
      double precision k, r0(6, nLp)
      complex*16 as(nLp, nmp), bs(nLp, nmp),as0(nLp, nmp),bs0(nLp, nmp),
     $           p0(nLp, nmp), q0(nLp, nmp),pp0(nLp, nmp),qq0(nLp, nmp),
     $           aMie(nLp,np), bMie(nLp,np),atr0(ni0,nij),btr0(ni0,nij), 
     $           A, ci, cin, cmz

      common/fnr/    fnr(0:2*(np+2))

      ci     = dcmplx(0.d0, 1.d0)
      cin    = dcmplx(0.d0,-1.d0)

C     should including the position of dipole, 
C     remember that nL = #particle + 1 when reading from dat file
      do i = 1,nL
        ind(i)=0
        do imn=1,nmp
            p0(i,imn)=0.d0
            q0(i,imn)=0.d0
            as(i,imn)=0.d0
            bs(i,imn)=0.d0

            as0(i,imn)=0.d0
            bs0(i,imn)=0.d0
         enddo
      enddo
C
c     This describes the field of a dipole oriented along the x-axis and
c     located at the position of the sphere nL in the basis { N3(nL), M3(nL) }
C
      as0(nL, 1) = ci/fnr(3)
      as0(nL, 3) = -ci/fnr(3)

C     展开dipole到as,bs中
c     Tranform the dipole field from the basis { N3(nL), M3(nL) } into
c     the bases { N1(i), M1(i) | i < nL }. The result is stored in as, bs
      call trans(nL, r0, nmax, uvmax, fint, atr0, btr0, ek, drot, as0,
     $     bs0, as, bs, ind)

c     multiply the as(m,n) and bs(m,n) by [- aMie(n)] and copy the result
c     to the p0(m,n) and q0(m,n)
c     We need to switch sign because the expansion of the incident
c     field has a leading minus, while that of the dipole field does not
c
      nL = nL - 1
      do i = 1,nL
         write(6, *) 'Sphere ', i
         do n=1,nmax(i)
            imn = n*n
            do m=-n,n
C      save the original incident dipole
                pp0(i,imn) = - as(i,imn) * 1.0D0
                qq0(i,imn) = - bs(i,imn) * 1.0D0
C     set the value for interaction solver
                p0(i,imn) = - as(i,imn) * aMie(i,n)
                q0(i,imn) = - bs(i,imn) * bMie(i,n)
C     initial guess
                as(i,imn) = p0(i,imn)
                bs(i,imn) = q0(i,imn)
c     write(6, *) n, m, as(i,imn), p0(i,imn)
c     write (6, *) i, imn, pp0(i,imn), qq0(i, imn)
                imn = imn + 1
            enddo
         enddo
      enddo

      return
      end
